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^ ■ Abstract 

o. 

'nJ" . Self-sustained subthreshold oscillations in a discrete-time model 

of neuronal behavior are considered. We discuss bifurcation scenar- 
io ■ ios explaining the birth of these oscillations and their transformation 
^ . into tonic spikes. Specific features of these transitions caused by the 

qh| discrete-time dynamics of the model and the influence of external noise 

L^ ' are discussed. 

1 Introduction 

Studies of dynamical behavior of biological networks require numerical sim- 
ulations of arrays containing a very large number of neurons. Despite the 
variety of physiological processes involved in the formation of neuron activity, 
the thorough studies of the large-scale networks need simple phenomenolog- 
ical models that can replicate the dynamics of individual neurons. Various 



suggestions for the design of low-dimensional maps for modeling the neu- 
rons' behavior have been proposed, see for example [H I2l El IH El El and 
references therein. Most of them were focused on the replication of either 
fast spikes or relatively slow bursts while the mechanisms for generation of 
specific footprints of spikes were neglected. 

A simple discrete-time model replicating the spiking-bursting neural ac- 
tivity has been suggested recently in [7|. This model is a 2D map that mimics 
rather realistically various types of transitions that occur in biological neu- 
rons. These transitions include routes between silence and tonic spiking as 
well as a triplet: silence <-h> bursts of spikes ^-^ tonic spiking. Such simple 
phenomenological models bear a high potential for further developments of 
computationally efficient methods for studies of functional behavior in large- 
scale neurobiological networks |H]. 

The bifurcation analysis of the map model carried out in |H] has shown 
that the transition from silence (a stable fixed point) to generation of ac- 
tion potentials is characterized by a sub-critical Andronov-Hopf bifurcation 
when an unstable invariant closed curve collapses into the stable fixed point. 
Therefore, the original map-model [7| provides only an abrupt transition 
from silence to spiking as a control parameter (e.g. the depolarization cur- 
rent) passes the excitability threshold. This scenario is quite typical for most 
types of biological neurons. However, experimental studies suggest that some 
neurons may come out of the silence softly through the regime of small oscil- 
lations below the threshold of the spike excitation ^0]. These subthreshold 
oscillations of almost sinusoidal form facilitate the generation of spike os- 
cillations when the membrane gets depolarized or hyperpolarized [TTl IT^ . 
These small oscillations can play an important role in shaping specific forms 
of rhythmic activity that are vulnerable to the noise in the network dynam- 
ics [niiini- 

In this paper we modify the map model so that it can generate stable 
subthreshold oscillations. We start the discussion of the model dynamics 
with the analysis of local bifurcations of a fixed point of the map. We show 
that the loss of stability of the fixed point is accompanied by the birth of the 
stable invariant circle which initiate a family of canards in the map. Further 
evolution of the circle leads to a breakdown of the invariant circle that gives 
rise to chaos. We also elaborate on the mechanism of the onset of irregular 
spiking which is due to heteroclinic-like crossings between the stable and 
unstable invariant sets, which are the images of the slow motion "surfaces" 
in the unperturbed map. The role played by small subthreshold oscillations 



in the responsiveness of the map to external noise is considered. 

2 Map-based model with stable subthreshold 
oscillations 

The map-based model of spiking-bursting neuron oscillations, following [7], 
can be written in the form of the two-dimensional map 

^ = faix,y + /3), (la) 

y = y- Kx + l-a), (lb) 

where the x-variable replicates the dynamics of the membrane potential, 
the parameters a, a and < // -C 1 control individual dynamics of the 
system. Some input parameters P and a are employed to provide coupling 
with other such models afterwards; both stand for injected currents. The 
principal distinction of the original map analyzed in [71 |H] and the one in 
question is camouflaged in the function fa{x,y + /3), which is given by 

{-a'^/4- a + y + f3, x < -1 - a/2, 

ax + {x + iy + y + P, -l-a/2<x<0, 

y + l+P, 0<x<y + l + (], ^^ 

-1, x>y + l + /3. 

The graph of this function is pictured in Fig. ^ The shape of nonlinear 
function is meant to achieve a replication of sharp tonic spikes in the dy- 
namics of the x-variable in the map. The slow ^/-variable can turn the spike 
generator on and off. The main difference between the function Q and that 
in the map proposed in [7] is its shape on the left hand side of the disconti- 
nuity point, i.e. a.tx<y + l + l3. Now the function contains an interval of 
parabola instead of a hyperbola used in [THH]. This modification is crucial for 
stability of subthreshold oscillations. The parabola reaches its minimum at 
X = —1 — a/2. At this point the graph of the function is continued leftward 
by a horizontal line, see Fig. ^ 

2.1 Fast and slow motions of the map 

For map (1), when // = 0, the slow subsystem (fTB]) is decoupled from the fast 
subsystem (llaj) . in which y is regarded as a perturbation parameter. One 




Figure 1: Geometry of the fast subsystem, map (fTaj) . for the parameter values 
at the tangent bifurcation: a = 1 and y = 0. The function is discontinuous 
at the point x = y + 1 + (3 which belongs to the rightmost interval, see (j21). 



may see from Fig. Qthat depending on y, the fast subsystem may have two 
fixed points, one stable and one unstable, or no fixed points. The transition 
between these states occurs via a saddle-node bifurcation. When y is varied, 
the fixed points trace out a parabola in the {y, x)-plane as shown in Fig. |21 
The traces of stable and unstable fixed points form on the {x,y) stable, Sps, 
and unstable, Spu, branches, respectively. 

The point of intersection of these branches with the nuUcline of slow 
subsystem flTB]) . which is given by Xfp = a — 1, is a fixed point of the two- 
dimensional map (1). It is easy to see that this fixed point O is stable if it 
is located on Sps and is unstable if it is on Spu- The case where the nuUcline 
crosses the parabola at the fold point requires a more delicate analysis. 

The nonstandard analysis ^7j predicts that when < /x <C 1, the nor- 
mally hyperbolic branches Sps and Sps persist in the form of the stable and 
unstable "slow motion" sets Ss and Su that remain /^-close to the originals, 
but yU^/^-close to them near the fold |18j . 

The idea behind generation of tonic spikes in the map (1) is illustrated 
in Fig. 121 The period of spiking is comprised of the following phases: the 
rest phase, where the phase point slides along the set Sg at a rate of order 
/i towards the fold point. Then, the phase point jumps up indicating the 
beginning of a spike. Its upward motion is stopped by a delimiter (the third 
segment of function Q) that reflects the phase point towards the stable slow 
surface Ss through the line segment x = —1. 
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Figure 2: (y, x) -bifurcation diagram of the fast subsystem, map (fTaj) . for 
a = 1 and /i = 0. The branches Spu and Sps are traced out by the stable 
and unstable fixed points of the fast map as the parameter y varies. When 
< yU -C 1, this {y,x) plane becomes the phase portrait of the 2D map. 
The shown trajectory of the 2D map is computed at /i = 0.04, a = 0.99 and 
a = 0. The dashed curves are the zero level ones for function (J2|. 



2.2 Birth of invariant curve 

Next we carry out the bifurcation analysis of the fixed point O. Our con- 
sideration is restricted to the domain a^/2 — a + y + j3<x<0, i.e. to 
the parabolic segment of function (J2)). This means that the fast subsys- 
tem (fTaj) is chosen to be close to the tangent bifurcation. The moment of 
the bifurcation is pictured in Fig. ^ From (jlbjl one finds the x-coordinate 
of the fixed point Xfp = a — 1. Therefore, the fixed point is located at 
the parabolic segment of function when the parameter values are within the 
range —a/2 < a < 1. The y-coordinate of the fixed point for this range of 



parameters is yjp = {a — 1)(1 



a 



p. 



Since fixed point O is a single fixed point of map (1), the further stability 
analysis is reduced to two possible local bifurcations: a flip where one of the 
multipliers of the fixed point equals —1; and the Andronov-Hopf bifurcation, 
where the fixed point possesses a pair of multipliers equal e''^*'^ on the unit 
circle. Simple calculations reveal that the flip bifurcation takes place outside 
of the considered parameter region. 

In the case of the Andronov-Hopf (AH) bifurcation, the Jacobian 

a + 2a 1 
— /i 1 



of the map equals 1 at the fixed point, while its trace equals 2cos(y9. The 
equation of the corresponding bifurcation curve AH is given by 

aAH = -2cr + 1 - /i. (3) 

On this curve, the fixed point has a pair of complex conjugate multipliers: 



Pi,2 = i - - ± -^/^i{4:- fi) = cosij±ismip. 



(4) 



To determine the stability of fixed point O right at the bifurcation state we 
need to evaluate the sign of the first Lyapunov coefficient Li. Note that the 
value of Li is not an invariant as it depends on coordinate transformations. 
The critical fixed point O is stable if Li < 0, and unstable if Li > 0. 
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Figure 3: Bifurcation diagram illustrating the birth of "small" subthreshold 
oscillations transforming into spikes as parameter a increases. Top and bot- 
tom branches corresponds to the highest and lowest values of the x-variable 
for a given value of parameter a. The other parameters of the map are set 
as follows a = 0.99, /3 = 0.0 and /i = 0.02. 

Let us introduce new coordinates in which the fixed point is translated 
to the origin: 



X + 1 — a 
y + (a - l){a - I) + a^ + p 



Now the map looks as follows 

a; = (1 — fi)x + x'^ + y, V = V — f^x. 



(5) 
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Let us next make another transformation 

x\ f 1 \U 

y J y sin ijj 1 — cos ip J \ rj 

with sinip and cosip defined in (0}. This makes the hnear part of (0) a 
rotation through ip: 

e ^ / COST/; -sin^A, W ^ V /^ ^1^ ^ 2^ 
7] J y sin Ip cosijj, J \ V J V -'- / 

In variable z = C, + if] the map assumes the complex form 



with 



# I "--iU 2 , * I '-'02 *2 /c\ 

26 ^ + —2; + ciizz + —z (6) 

Zi Zi 



C20 = 2 ( ^^'^ ^ ~ ' ^11 " 2 ( ^ ~ ^^"^ 2" ) ' ^02 = 2 ( ta'^ Y - • C^) 
The normalizing transformation 

C20 2 '^11 * C02 *2 

~ ■ ■ ~ -Z — —ZZ r-- --Z 



eliminates all the quadratic terms in (jH)) so that the coefficient Li at the 
desired cubic term in the resulting normal form 

z = e'^z{l + (Li + tSi)zz*) + Odl^in 

becomes the sought Lyapunov value. Its expression reads as follows: 

_ (l-2e.*)e-^^c.„c„ _ ^ _ K^P, («, 

2(1 - e'^) 2 4 ^ ^ 

Substituting (|7j) into (jH)) yields 

^ 1 cos(^) _ 2 - /i ^ ^ 

^ 4 cos(^) + 1 4(4 - /i) 

for all small /i. Thus, the loss of stability of fixed point O on the Andronov- 
Hopf bifurcation curve is accompanied by the birth of a stable closed invariant 
curve emerging from O. This mechanism of the birth of subthreshold periodic 
oscillations is illustrated in Fig. |S1 as the parameter a increases. 
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Figure 4: Four waveforms of iterates x„ versus the discrete time n show 
the influence of an external Gaussian white noise with standard deviation s 
on the map oscillations. The control parameters of the map are a = 0.99, 
a = -0.0001 and /i = 0.02. 

2.3 Noise and subthreshold oscillations 

It was observed recently that subthreshold periodic activity shapes stochas- 
tic properties of spiking in a neuron influenced by noise, see for exam- 
ple O CH 113 Cni UH 120] • To illustrate similar properties in our model 
let us now consider the regime of subthreshold oscillations in the map with 
the presence of noise. For the sake of briefness we consider only the case when 
noise is applied to the fast subsystem. The discussion of the effects of noise, 
occurring in the fast and slow subsystems, can be found elsewhere J2I1 122] • 
The system (1) in presence of noise can be written as follows 



yn+1 = i/„ -/i(x„ + 1 - a). 



(9) 
(10) 



where (n is a delta-correlated Gaussian White Noise (GWN) with zero mean 
value and standard deviation value s. 

Waveforms of map model (1) operating in the regime of subthreshold 
oscillations and the influence of noise are shown in Fig. |3] Four traces plotted 
in the figure correspond to different values of the standard deviation, s, of 
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ISI [n] 



Figure 5: Histograms of the interspike interval (ISI) distribution computed 
for the waveforms shown in Fig. |3J (a) - s = 0.0002; (b) - s = 0.002 and (c) 
- s = 0.02. 



GWN. The top trace presents the case where the level of noise is insufficient 
to induce an action potential (a spike). When the level of noise exceeds a 
critical level the map starts producing occasional spikes. Such spikes are 
more likely to occur at the top of oscillation. This behavior results in the 
formation of a multi-hump structure in the probability distribution function 
of interspike intervals (ISIs), see histograms shown in Fig. |Sfa,b). 

A further raise of the noise level increases the probability of the action 
potentials. As the result, spikes occur almost every period of the subthreshold 
oscillation (see the bottom trace in Fig. |3)) and the distribution of probability 
density of ISIs transforms into a single hump structure, see Fig. |3fc). These 
results are in good agreement with the results obtained from ODE based 
neuron models [121 120] ■ 



3 Tangles of critical curves and chaos 



In the numerical simulations of map (1) we found that subthreshold oscilla- 
tions may be interrupted by irregular spiking even in the absence of external 
noise. An example of such a behavior is presented in Fig. IHl This interme- 
diate dynamics is observed only within a rather thin parameter interval at 
the border between the regimes of continuous subthreshold oscillations and 
tonic spike generation. To understand the dynamical mechanisms behind 
this sporadic spiking we numerically studied the evolution of an invariant 
circle as the parameter values enter this thin region. 
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Figure 6: Sporadic spiking in the map (1) computed with a = 1.25, a = 
—0.13 and fi = 0.02. The phase portrait of these oscillations is shown in 
Fig. Ufa). 

It follows from the theory of canards that the parameter domain for the 
existence of a stable invariant circle is a narrow strip of the order of 0{n) 
that adjoins the bifurcation curve AH. Furthermore, the size of the circle 
increases abnormally fast as the parameter values deviate from the Andronov- 
Hopf bifurcation curve and approach the critical values where the invariant 
circle breaks down. Such extreme sensitivity to slight parameter deviations 
hamper the detailed analysis of bifurcations associated with the circle as it 
breaks. 

The breakdown of the stable invariant circle leads to an interesting sit- 
uation depicted in Fig. [7|^a). One can observe the co-existence of two kinds 
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Figure 7: (a) Chaos in the noiseless map (1) computed with a = 1.25, 
a = —0.13 and fi = 0.02 The corresponding waveform is shown in fig. [Ul 
(b) Forward iterates of a small interval of the stable critical set Sg reveal 
increasing wiggles occurred around the unstable critical set Su- 

of special solutions following the unstable slow branch S^- They are called 
canards with a head (a spike) and ones without it J7j. A canard is charac- 
terized by a growing level of exponential instability with respect to nearby 
solutions. This instability is a necessary first component of chaotic behav- 
ior observed in a system. In addition, the presence of two types of canards 
creates mixing and uncertainty, which is the second important ingredient for 
the onset of chaos. This type of chaos in the neuron model (1) appears as 
small subthreshold oscillations alternating with sporadic spikes, see Fig. |B] 

To understand the dynamical mechanism behind the splitting of canards 
into two types, we conducted a numerical analysis of the behavior of Sg 
and Su for various parameter values selected close to the threshold for the 
breakdown of the stable invariant circle. To plot Ss, we iterated forward a 
large number of phase points initiated from a relatively short interval on the 
stable branch Sg, chosen as an initial approximation. Figure [7|^b) shows how 
the connected forward images of this interval become non-smooth, generating 
growing wiggles around the unstable set Su- Note that the unstable set Sg 
is easily computed using inverse mapping , see Appendix A. The insert in 
Fig. U][h) shows the zoomed-in wiggles. One can see from Fig. [7| (a) and (b) 
that upon entering the wiggling area, the phase point can land back in the 
stable set Sg, thereby completing another round of subthreshold oscillations, 
or jump up to make a spike. Such a situation is referred to as dynamical 
uncertainty [221 • Loosely speaking, the place on Su where the wiggles become 
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noticeable indicates the threshold beyond which the behavior of the map 
becomes uncontrollable. 

This situation is similar to the case of a periodically driven pendulum 
near a homoclinic orbit associated with the saddle point at the top. Such 
a system is typically studied using a Poincare return mapping defined over 
the period of the external force. Under proper conditions, the saddle fixed 
point of the mapping will possess transversal crossings between its stable and 
unstable sets (formerly, the stable and unstable separatrices of the saddle 
equilibrium in the autonomous system). These crossings generate the Smale 
horseshoe and, therefore, symbolic shift-dynamics in the system, see details 
in [T71 I2SI 121] and references therein. By construction, this type of Poincaree 
mapping is a diffeomorphism. However, our map (j21) is an endomorphism, 
e.g. a non- invert ible one, and as a result it may possess some exotic features 
prohibited in the invertible maps, see more in [221 • O^^ ^^ those features is 
that the set Sg may self-cross (so does the set Su in the backward time). This 
situation is sketched in Fig. IHl^a) and also can be seen on the trajectories' 
behavior in the lower left corner of the attractor shown in Fig. [7||^a). One can 
assume that these self-crossings can stimulate the conditions for the onset of 
a topological Small horseshoe (see Fig. IHfb)), whose presence is a de-facto 
proof of complex chaotic dynamics. 
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Figure 8: Heteroclinic-like crossings of the critical sets 5* and S"". Self- 
crossings of the set S^, which are one of the features of noninvertible maps, 
may generate a topological Smale horseshoe. 
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4 Conclusion 

It is shown that a simple map can be employed to replicate the behavior of 
neurons with self-sustained subthreshold oscillations. These oscillations are 
achieved by a special selection of a nonlinear function in the fast subsystem. 
The dynamical mechanisms behind the transitions from silence to subthresh- 
old oscillations of small amplitude and then to spiking activity are explained 
using the bifurcation approach. 

Here we focused mostly on the individual dynamics of the map-based 
model. As a result we considered only the case when (3 and a are constants. 
One may notice from (1) that the parameter (3 can be eliminated by using 
the variable transformation y + f] ^—>- y. However, the role of input parameter 
/3 becomes important when a time dependant input is considered, see j7] for 
detail. We would like to note that for studies of non-autonomous dynamics of 
this map model one needs to modify function Q to insure that no trajectory 
of (1) gets locked in the interval 0<x<y + l+(3a.s(3oTy increase. The 
discussion on suggested alterations of the function to resolve this problem 
can be found elsewhere jTj. 
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6 Appendix: Inverse map 

To locate the unstable set iS" which is a surface of slow motion one should 
consider the inverse map Finv '■ {x,y) \-^ {x,y) G D defined in D := { — 1 — 
a/2 < X < 0}. Within D the inverse F^v assumes the following form 

X = ax + {x + lf + y + P, (11a) 

y = y- jJ-ix + l-a). (lib) 

Subtracting ()llb|l from (jlla|) one gets 

x-y = ax + {x + l)"^ + H{x + l-a) + p. (12) 

Solving it for x one can find x as a the following function in x and y 

x = -{2 + a + /i)/2 + ^{2 + a + /i)74 - {1 + (3 + fi{l - a) - x + y. (13) 

To get the equation for the ^/-variable expression (fT!?|l must be plugged into 
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